#################################################################################
# Replication file
# Where to flee? Preferences for host communities among displaced people in Congo
# For questions: sigrid.weber@ie.edu
# Last updated: September 24, 2025 - SW
#################################################################################

# setup ------------------------------------------------------------------------

# load required packages
library(tidyverse) # for data handling
library(estimatr) # for robust SEs and other models
library(cregg) # for conjoint analyses 
library(here) # for easier/relative/automatic filepaths 
library(patchwork) # assembling plots
library(modelsummary) # nice table outputs for regressions
library(legendry) # nested labels
library(corrplot) # correlation matrix
library(vtable)
library(likert)
library(cjoint)

# set file path to directory of the replication R code
here("<your file path>")

# read data in relative file path

# conjoint data
conjoint <- read_csv(here("replication_data/where2flee_conjoint.csv"))

# full study sample 
all <- read_csv(here("replication_data/where2flee_data.csv"))

# sample that received conjoint
data <- all|> filter(ID_Household %in% conjoint$ID_Household)


# prepare functions for easy plotting ------------------------------------------

# write function to add attributes to the plot and enforce order of estimates
beautify_plot <- function(data){
  # how often do labels have to repeat?
  num_rep = if_else(is.null(data$BY),1,length(unique(data$BY)))
  
  # make nice labels for features
  data = data |> 
    mutate(feature = case_when(feature == "attr_economic" ~ "Economic \nintegration",
                               feature == "attr_social" ~ "Social \nintegration",
                               feature == "attr_local" ~ "Relational \nintegration",
                               feature == "attr_cultural" ~ "Cultural \nintegration",
                               feature == "attr_political" ~ "Political \nintegration"),
           # order the features in the same way
           order = rep(rep(5:1, each = 2), num_rep*2))
  
  data
}

# write a function to plot the different models 
plot_conjoints <- function(input, estimand= "marginal mean",subgroup,plot_type="side"){
  
  clean_input = beautify_plot(input)
  int = if_else(estimand == "marginal mean",0.5,0)
  
  if(is.null(clean_input$BY)){
    ggplot(clean_input, 
           aes(estimate,  reorder(paste0("   ",level, "&", feature),order))) +
      geom_point(size = 2) + 
      guides(y = guide_axis_nested(key=key_range_auto(sep="&"))) +
      theme_bw() +
      geom_errorbarh(aes(xmin = lower, xmax = upper, height=.01,
                         linewidth= alpha)) + 
      scale_linewidth_continuous(range = c(0.5,1.1))+
      geom_vline(xintercept = int, linetype = "dashed", color= "gray60") +
      labs(x = "", y = "") + 
          theme(text=element_text(size=13),
            plot.title = element_text(size=14))+
      guides(linewidth= "none")
    
    
  } else{
    
    if(plot_type == "side"){
      ggplot(clean_input, 
             aes(estimate,  reorder(paste0("   ",level, "&", feature),order))) +
        geom_point(size = 2) + 
        guides(y = guide_axis_nested(key=key_range_auto(sep="&"))) +
        theme_bw() +
        geom_errorbarh(aes(xmin = lower, xmax = upper, height=.01,
                           linewidth= alpha)) + 
        scale_linewidth_continuous(range = c(0.5,1.1))+
        geom_vline(xintercept = int, linetype = "dashed", color = "gray60") +
        labs(x = "", y = "") + 
        facet_wrap(~ BY, ncol = 2)+
        theme(text=element_text(size=13),
              plot.title = element_text(size=14))+
        guides(linewidth = "none")
    } else{
      
      ggplot(clean_input, 
             aes(estimate,  reorder(paste0("   ",level, "&", feature),order),color=BY)) +
        geom_point(size = 2,position = position_dodge(width = .5)) + 
        guides(y = guide_axis_nested(key=key_range_auto(sep="&"))) +
        scale_color_grey()+
        theme_bw() +
        geom_errorbarh(aes(xmin = lower, xmax = upper, height=.01,
                           linewidth= alpha),
                       position = position_dodge(width = .5)) + 
        scale_linewidth_continuous(range = c(0.5,1.1))+
        geom_vline(xintercept = int, linetype = "dashed", color = "gray60") +
        labs(x = "", y = "", color = str_to_sentence(subgroup)) + 
        theme(text=element_text(size=13),
              plot.title = element_text(size=14,hjust = 0.5),
              legend.position = "bottom")+
        guides(linewidth = "none")
      
      
    }
  }
  
}

# descriptives -----------------------------------------------------------------

# number of households in conjoint
length(unique(conjoint$ID_Household))

# demographics of main respondent
round(100*prop.table(table(data$pop_group, useNA = "always")),2)
round(100*prop.table(table(data$gender, useNA = "always")),2)
round(100*prop.table(table(data$married, useNA = "always")),2)
round(100*prop.table(table(data$schooled, useNA = "always")),2)
round(100*prop.table(table(data$graduated, useNA = "always")),2)
round(100*prop.table(table(data$literacy, useNA = "always")),2)
round(100*prop.table(table(data$tshiluba,useNA = "always")),2)
round(100*prop.table(table(data$working, useNA = "always")),2)
round(100*prop.table(table(data$regular_income, useNA = "always")),2)
round(100*prop.table(table(data$community_meeting, useNA = "always")),2)
round(100*prop.table(table(data$christian, useNA = "always")),2)
round(100*prop.table(table(data$religion, useNA = "always")),2)

# main effect ----------------------------------------------------------

# set right factor reference
conjoint <- conjoint |> 
  mutate(across(starts_with("attr_"),factor)) |> 
  mutate(attr_economic = relevel(attr_economic, ref = "Difficult to get work in community"),
         attr_social = relevel(attr_social, ref = "Cannot join local church"),
         attr_local = relevel(attr_local, ref = "No relatives in the community"),
         attr_cultural = relevel(attr_cultural, ref = "Few speak mother tongue"),
         attr_political = relevel(attr_political, ref = "Not allowed at village meetings"))

# amce model
main_model_amce <- bind_rows(cj(chosen ~ 
                                  attr_economic + attr_social + attr_local + 
                                  attr_cultural + attr_political, 
                                id = ~ID_Household, data = conjoint,
                                estimate = "amce", alpha = 0.05) |> 
                               mutate(alpha=0.05),
                             cj(chosen ~ 
                                  attr_economic + attr_social + attr_local + 
                                  attr_cultural + attr_political, 
                                id = ~ID_Household, data = conjoint,
                                estimate = "amce", alpha = 0.1) |> 
                               mutate(alpha=0.1))

# results table
knitr::kable(main_model_amce |> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) , 
             format="latex",booktabs=T, digits=3)

# plot the main effects 
plot_conjoints(main_model_amce, estimand = "AMCE")

# uncomment to save plot in figures folder
# ggsave("figures/main_amce_plot.jpg", width= 8, height= 5)


# formal test of differences
model_for_test <- cj(chosen ~ 
                       attr_economic + attr_social + attr_local + 
                       attr_cultural + attr_political, 
                     id = ~ID_Household, data = conjoint,
                     estimate = "amce", alpha = 0.05) |> 
  filter(estimate != 0)


pairwise_test <- function(x1,x2){
  
  beta1 = model_for_test$estimate[x1]
  beta2 = model_for_test$estimate[x2]
  se1 = model_for_test$std.error[x1]
  se2 = model_for_test$std.error[x2]
  
  # calculate differences
  diff <- beta1 - beta2
  se_diff <- sqrt(se1^2 + se2^2)
  
  # t statistic
  t_val = diff/se_diff
  
  # p value
  p_val = 2 * (1 - pnorm(abs(t_val)))
  
  p_val
}


formal_test <- data.frame(
  z = c("Econ","Social","Local","Cultural","Political"),
  Econ = c(pairwise_test(1,1), pairwise_test(2,1), pairwise_test(3,1),pairwise_test(4,1),pairwise_test(5,1)),
  Social = c(pairwise_test(1,2), pairwise_test(2,2), pairwise_test(3,2),pairwise_test(4,2),pairwise_test(5,2)),
  Local = c(pairwise_test(1,3), pairwise_test(2,3), pairwise_test(3,3),pairwise_test(4,3),pairwise_test(5,3)),
  Cultural = c(pairwise_test(1,4), pairwise_test(2,4), pairwise_test(3,4),pairwise_test(4,4),pairwise_test(5,4)),
  Political = c(pairwise_test(1,5), pairwise_test(2,5), pairwise_test(3,5),pairwise_test(4,5),pairwise_test(5,5))
) |> 
  mutate(across(-z,\(x)if_else(x==1,NA,round(x,3)))) |> 
  mutate(across(-z, \(x) if_else(x<0.05,paste0(x,"*"), as.character(x)))) 

kable(formal_test |> rename(`X`= z), 
             format="latex",booktabs=T)


# mechanisms ------------------------------------------

ranking_econ_contribution <- beautify_plot(bind_rows(
  cj(ranking_economic_contribution ~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.05) |> mutate(alpha=0.05),
  cj(ranking_economic_contribution ~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.1)|> mutate(alpha=0.1)))
  
ranking_feeling_welcome <- beautify_plot(bind_rows(
  cj(ranking_feeling_welcome~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.05) |> mutate(alpha=0.05),
  cj(ranking_feeling_welcome ~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.1)|> mutate(alpha=0.1)))

ranking_feeling_safe<- beautify_plot(bind_rows(
  cj(ranking_feeling_safe~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.05) |> mutate(alpha=0.05),
  cj(ranking_feeling_safe~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.1)|> mutate(alpha=0.1)))

ranking_contribute_ideas <- beautify_plot(bind_rows(
  cj(ranking_contribute_ideas~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.05) |> mutate(alpha=0.05),
  cj(ranking_contribute_ideas~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.1)|> mutate(alpha=0.1)))

ranking_trust_community <- beautify_plot(bind_rows(
  cj(ranking_trust_community~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.05) |> mutate(alpha=0.05),
  cj(ranking_trust_community~ attr_economic + attr_social + attr_local + 
       attr_cultural + attr_political, 
     id = ~ID_Household, data = conjoint,estimate = "amce",alpha=0.1)|> mutate(alpha=0.1)))

ranking_mechanisms <- bind_rows(ranking_econ_contribution,
                                ranking_feeling_welcome,
                                ranking_feeling_safe,
                                ranking_contribute_ideas,
                                ranking_trust_community) |> 
  mutate(outcome = case_when(
    estimate == 0 ~ "Reference category",
    outcome == "ranking_economic_contribution" ~ "Economic contribution",
    outcome == 'ranking_feeling_welcome' ~ "Feeling welcome",
    outcome == "ranking_feeling_safe" ~ "Feeling safe",
    outcome == "ranking_contribute_ideas" ~ "Contribute ideas",
    outcome == "ranking_trust_community" ~ "Trust community"
  )) |> 
  mutate(outcome = relevel(factor(outcome),ref="Reference category"))

ggplot(ranking_mechanisms,
       aes(estimate,  reorder(paste0("   ",level, "&", feature),order),
           color=outcome, shape=outcome)) +
  geom_point(size = 3,position = position_dodge(width = .7)) + 
  guides(y = guide_axis_nested(key=key_range_auto(sep="&"))) +
  theme_bw() +
  geom_errorbarh(aes(xmin = lower, xmax = upper, height=.01,
                     linewidth = alpha),position = position_dodge(width = .7)) + 
  scale_linewidth_continuous(range = c(0.5,1.1),guide="none")+
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
  labs(x = "", y = "", color="", shape="") + 
  scale_color_grey()+
  theme(text=element_text(size=13),
        legend.position = "bottom",
        plot.title = element_text(size=10))+
  guides(color=guide_legend(ncol=4),shape=guide_legend(ncol=3))

knitr::kable(ranking_mechanisms|> 
    select(-statistic, -feature,-order) |> 
    filter(is.na(std.error)==F)|> 
      pivot_wider(id_cols=c(outcome,level,estimate,std.error,z,p),
                  names_from = alpha,
                  values_from = c(lower,upper)), 
  format="latex",booktabs=T, digits=3)

# uncomment to save plot in figures folder
# ggsave("figures/amce_by_mechanisms.jpg", width=8, height= 5)

# subgroup analysis ----------

conjoint <- left_join(conjoint, data)

# Displacement status 2 - based on pop_group2 (i.e., "host community" vs "others")
conjoint <- conjoint |> 
  mutate(pop_group2 = factor(pop_group2),
         pop_group = factor(pop_group),
         tshiluba_speaker = factor(tshiluba, levels = c(0,1), labels= c("Other mother tongue","Tshiluba speaker")),
         can_read=factor(literacy))

# by all 4 subgroups
subgroupx <- bind_rows(cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint,
                          by= ~pop_group, 
                          estimate = "amce", alpha = 0.05) |> 
                         mutate(alpha=0.05),
                       cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint, by= ~pop_group,
                          estimate = "amce", alpha = 0.1) |> 
                         mutate(alpha=0.1))

plot_conjoints(subgroupx, estimand = "amce", 
               subgroup = "Displacement groups",plot_type = "color")+
  guides(color = guide_legend(nrow = 2))
#ggsave("figures/subgroup_displacement_status.jpg", width= 8.5, height= 5)

# results table
knitr::kable(subgroupx |> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(BY, level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) , 
             format="latex",booktabs=T, digits=3)


# test of difference between subgroups
amce_diffs(conjoint, chosen ~ attr_economic + attr_social + attr_local + attr_cultural + attr_political, 
   id = ~ID_Household, by= ~pop_group) |> filter(p<0.05)|> 
  select(-outcome,-statistic, -feature,-pop_group) |> 
  mutate(across(estimate:upper, \(x) round(x,3))) |> 
  kable(format ="latex",booktabs=T)


# appendix figures: additional subgroups -------

# by displaced or not 
subgroupx <- bind_rows(cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint,
                          by= ~pop_group2, 
                          estimate = "amce", alpha = 0.05) |> 
                         mutate(alpha=0.05),
                       cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint, by= ~pop_group2,
                          estimate = "amce", alpha = 0.1) |> 
                         mutate(alpha=0.1))

plot_conjoints(subgroupx, estimand = "amce", subgroup = "",plot_type = "color")
#ggsave("figures/subgroup_displacement_status2.jpg", width= 8, height= 5)

# results table
knitr::kable(subgroupx |> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(BY, level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) , 
             format="latex",booktabs=T, digits=3)

# Tshiluba speaker
subgroupx <- bind_rows(cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint,
                          by= ~ tshiluba_speaker, 
                          estimate = "amce", alpha = 0.05) |> 
                         mutate(alpha=0.05),
                       cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint, by= ~ tshiluba_speaker,
                          estimate = "amce", alpha = 0.1) |> 
                         mutate(alpha=0.1))

plot_conjoints(subgroupx, estimand = "amce", 
               subgroup = "Mother tongue",plot_type = "color")
#ggsave("figures/subgroup_tshiluba.jpg", width= 8.5, height= 5)


# Literacy
subgroupx <- bind_rows(cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint,
                          by= ~ can_read, 
                          estimate = "amce", alpha = 0.05) |> 
                         mutate(alpha=0.05),
                       cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint, by= ~ can_read,
                          estimate = "amce", alpha = 0.1) |> 
                         mutate(alpha=0.1))

plot_conjoints(subgroupx, estimand = "amce", 
               subgroup = "Literacy",plot_type = "color")
#ggsave("figures/subgroup_literacy.jpg", width= 8.5, height= 5)


# appendix figures: marginal means ---------------------------------------------

# marginal means overall
main_model_mm <- bind_rows(cj(chosen ~ 
                                attr_economic + attr_social + attr_local + 
                                attr_cultural + attr_political, 
                              id = ~ID_Household, data = conjoint,
                              estimate = "mm", alpha = 0.05) |> 
                             mutate(alpha=0.05),
                           cj(chosen ~ 
                                attr_economic + attr_social + attr_local + 
                                attr_cultural + attr_political, 
                              id = ~ID_Household, data = conjoint,
                              estimate = "mm", alpha = 0.1) |> 
                             mutate(alpha=0.1))

plot_conjoints(main_model_mm, estimand = "marginal mean")

# results table
knitr::kable(main_model_mm |> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) |> 
               select(-p), 
             format="latex",booktabs=T, digits=3)


# uncomment to save plot in figures folder
# ggsave("figures/main_mm_plot.jpg", width= 8, height= 5)


# Marginal means for displacement status
subgroupx <- bind_rows(cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint,
                          by= ~pop_group, 
                          estimate = "mm", alpha = 0.05) |> 
                         mutate(alpha=0.05),
                       cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political, 
                          id = ~ID_Household, data = conjoint, by= ~pop_group,
                          estimate = "mm", alpha = 0.1) |> 
                         mutate(alpha=0.1))

plot_conjoints(subgroupx, estimand = "marginal mean", subgroup = "Displacement status",
               plot_type = "color")
# ggsave("figures/subgroup_displacement_status_mm.jpg", width=10.5, height= 5)

knitr::kable(subgroupx |> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(pop_group, level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) |> 
               select(-p), 
             format="latex",booktabs=T, digits=3)


# appendix figures: continuous outcome ---------------------------------------------

model_cont <- bind_rows(cj(ranking ~ 
                                attr_economic + attr_social + attr_local + 
                                attr_cultural + attr_political, 
                              id = ~ID_Household, data = conjoint,
                              estimate = "amce", alpha = 0.05) |> 
                             mutate(alpha=0.05),
                           cj(ranking ~ 
                                attr_economic + attr_social + attr_local + 
                                attr_cultural + attr_political, 
                              id = ~ID_Household, data = conjoint,
                              estimate = "amce", alpha = 0.1) |> 
                             mutate(alpha=0.1))

plot_conjoints(model_cont, estimand = "AMCE")

knitr::kable(model_cont|> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) , 
             format="latex",booktabs=T, digits=3)

# uncomment to save plot in figures folder
# ggsave("figures/main_cont.jpg", width= 8, height= 5)

# summary statistics -----------------------------------------------------------

sum_stats <- data |> select(ID_Household, pop_group,food_insecurity, food, ptsd,
                            trust_num, experience_violence, concern_violence,
                            malaria, married, gender,regular_income, current_hosts, 
                            community_meeting,schooled, graduated, literacy, working,
                            tshiluba, bubindi, tshikete, lingala, 
                            christian,protestant, catholic, revivalist, apostolic,
                            kimbanguist, animist, muslim, religosity,
                            ) |> 
  mutate(across(!c(ID_Household,pop_group, food_insecurity,trust_num, regular_income, current_hosts, 
                   community_meeting,tshiluba, bubindi, tshikete, lingala, 
                   christian,protestant, catholic, revivalist, apostolic,
                   kimbanguist, animist, muslim),
                \(x) if_else(x %in%c("Above median food insecurity", 
                                     "Above median levels of PTSD",
                                     "High trust","Above median levels of trust",
                                     "Has experienced violence","Has malaria",
                                     "High concerns about violence","Has returned",
                                     "Displaced","Currently displaced","Male",
                                     "Married", "Has been to school",
                                     "Higher than primary school",
                                     "Can write or read", "Working","Once a week or more",
                                     "Animist","Don't know","Jehovah's Witnesses", "Muslim"),1,0)))

names(sum_stats) <- c("ID","pop_group","Food insecurity (numeric)", "Food insecurity (binary)",
                      "PTSD", "Trust across groups (numeric)",
                      "Experienced violence","Concerned about violence","Positive malaria test",
                      "Married","Male","Regular income","Currently hosting","Community meeting attendance",
                      "Attended school", "Ever graduated","Literacy","Working",
                      "Tshiluba speaker", "Bubindi speaker","Tshikete speaker","Lingala speaker",
                      "Any Christian denomination","Protestant","Catholic","Revivalist","Apostolic",
                      "Kimbanguist","Animist","Muslim",
                      "Weekly church attendance")

sumtable(sum_stats,
         summ=c('notNA(x)',
                'median(x)',
                'mean(x)',
                'sd(x)',
                'min(x)',
                'max(x)'),
         group="pop_group",
         summ.names=c("Obs.","Median","Mean","St.Dev.","Min","Max"),
         digits = 3, fit.page = '\\textwidth', 
         out = 'latex'
         )


custom_f <- function(x, y){
  result = aov(sum_stats |> pull(!!x)~pop_group, sum_stats)
  p = summary(result)[[1]][["Pr(>F)"]][1]
  
  data.frame("Variable" = x,"F test" = p,
             "Significant" = if_else(p<0.05,"*",""))
}

apply_to <- c(names(sum_stats |> select(-ID,-pop_group)))
map_dfr(apply_to,custom_f) |> 
  mutate(F.test =round(F.test,3))


sumtable(sum_stats |> 
           mutate(host = as.numeric(pop_group =="Host community"),
                  idp = as.numeric(pop_group == "IDP"),
                  repatriated = as.numeric(pop_group== "Repatriated refugee"),
                  returned = as.numeric(pop_group=="Returned IDP"),.after=1) |> 
           select(-pop_group),
         summ=c('notNA(x)',
                'median(x)',
                'mean(x)',
                'sd(x)',
                'min(x)',
                'max(x)'),
         summ.names=c("Obs.","Median","Mean","St.Dev.","Min","Max"),
         digits = 3, fit.page = '\\textwidth', 
         out = 'latex'
)



# appendix: selection into conjoint ----
selection <- list(
  "Conjoint assigned" = lm(conjoint_received ~  pop_group + food_insecurity+trust_num+
                             malaria+married+gender+schooled + regular_income+ ptsd+literacy+working, all)
)

modelsummary(selection, stars=c("*"=0.05),
            # "latex",
             booktabs=T)


# appendix: adjusting for covariates -----
conjoint <- conjoint %>% 
  mutate(regular_income = if_else(regular_income==1 , "Regular income","No regular income")) %>% 
  mutate(across(
    c( food, concern_violence,malaria, married,gender,
       schooled, regular_income, graduated, literacy, working, religion,
       tshiluba_speaker),
    \(x) factor(x)
  ))


main_model_amce_cov <- bind_rows(cj(chosen ~ 
                                  attr_economic + attr_social + attr_local + 
                                  attr_cultural + attr_political+ food + concern_violence+ 
                                  tshiluba_speaker+religion+malaria+ married+gender+schooled+
                                  regular_income+graduated+literacy+working, 
                                id = ~ID_Household, data = conjoint,
                                estimate = "amce", alpha = 0.05) |> 
                               mutate(alpha=0.05),
                             cj(chosen ~ 
                                  attr_economic + attr_social + attr_local + 
                                  attr_cultural + attr_political+ food + concern_violence+ 
                                  tshiluba_speaker+religion+malaria+ married+gender+schooled+
                                  regular_income+graduated+literacy+working, 
                                id = ~ID_Household, data = conjoint,
                                estimate = "amce", alpha = 0.1) |> 
                               mutate(alpha=0.1))

# results table
knitr::kable(main_model_amce_cov |> 
               select(-outcome,-statistic, -feature) |> 
               filter(is.na(std.error)==F)|> 
               pivot_wider(id_cols=c(level,estimate,std.error,z,p),
                           names_from = alpha,
                           values_from = c(lower,upper)) , 
             format="latex",booktabs=T, digits=3)

# plot the main effects 
plot_conjoints(main_model_amce_cov|> filter(grepl("attr_",feature)), estimand = "AMCE")
#ggsave("figures/main_effect_covariate_adj.jpg", width= 8, height= 5)


subgroupy <- bind_rows(cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political + food + concern_violence+ tshiluba_speaker+religion+
                            malaria+ married+gender+schooled+regular_income+graduated+literacy+working, 
                          id = ~ID_Household, data = conjoint,
                          by= ~pop_group, 
                          estimate = "amce", alpha = 0.05) |> 
                         mutate(alpha=0.05),
                       cj(chosen ~ 
                            attr_economic + attr_social + attr_local + 
                            attr_cultural + attr_political + food + concern_violence+ malaria+ married+gender+schooled+regular_income+graduated+literacy+working, 
                          id = ~ID_Household, data = conjoint, by= ~pop_group,
                          estimate = "amce", alpha = 0.1) |> 
                         mutate(alpha=0.1))

plot_conjoints(subgroupy %>% filter(grepl("attr_", feature)), 
               estimand = "amce", subgroup = "",plot_type = "color")+
  guides(color=guide_legend(nrow=2,byrow=TRUE))
#ggsave("figures/subgroup_covariate_adj.jpg", width= 8, height= 6)

# appendix: interactions between features -----
interaction_results <- cjoint::amce(
  chosen ~ 
    attr_economic + attr_social + attr_local + 
    attr_cultural + attr_political + 
    
    attr_economic:attr_social+
    attr_economic:attr_local+
    attr_economic:attr_cultural+
    attr_economic:attr_political+
    
    attr_social:attr_local+
    attr_social:attr_cultural+
    attr_social:attr_political+
    
    attr_local:attr_cultural+
    attr_local:attr_political+
  
    attr_cultural:attr_political,
  data=conjoint, cluster=TRUE, respondent.id="ID_Household")

# Print summary

plot(interaction_results,plot.theme=theme_bw()+theme(legend.position = "none"))
#ggsave("figures/conjoint_interactions.jpg", width= 7, height= 7)
